library(tidyverse)
argv=commandArgs()
ref = read_tsv(argv[6])
alt = read_tsv(argv[7])

autoLength=function(x,y){
  x1 = as.integer(str_split(x,',',simplify = TRUE)[1,][1])
  x2 = as.integer(tail(str_split(x,",",simplify = TRUE)[1,],1))
  dis = x2 - x1 - 1
  flank = as.integer((y - dis) / 2)
  if((y - dis) %% 2){
    left = x1 - flank + 1
    right = x2 + flank}
  else{
    left = x1 - flank + 1 
    right = x2 + flank -1}
  return(c(left-1,right))
}

if(nrow(ref)>0){
ref = ref %>% separate(col =sequence_name , into = c('chr','pos','MNV_id','ref','alt','MNV_rsid'), sep = "_") %>%  
  mutate(pos0=sapply(pos,function(x){return(autoLength(x,60)[1])},simplify = TRUE)) %>% 
  mutate(start=pos0+start,stop=pos0+stop) %>% 
  unite("tag",pos,start,stop,sep = '-') %>% 
  mutate(tag1=sapply(tag,function(x){
    p0 = str_split(x,"-",simplify = TRUE)[1,][1]
    p = str_split(p0,",",simplify = TRUE)[1,]
    start = str_split(x,"-",simplify = TRUE)[1,][2]
    stop = str_split(x,"-",simplify = TRUE)[1,][3]
    for(i in p){if(i>=start & i<=stop) return(1)}
    return(0)
  },simplify = TRUE)) %>% filter(tag1==1) %>% separate(col=tag,into=c('pos','start','stop'),sep = '-') %>% 
  unite("sequence_name",'chr','pos','MNV_id','ref','alt','MNV_rsid',sep = '_')
}

if(nrow(alt)>0){
alt = alt %>% separate(col =sequence_name , into = c('chr','pos','MNV_id','ref','alt','MNV_rsid'), sep = "_") %>%  
  mutate(pos0=sapply(pos,function(x){return(autoLength(x,60)[1])},simplify = TRUE)) %>% 
  mutate(start=pos0+start,stop=pos0+stop) %>% 
  unite("tag",pos,start,stop,sep = '-') %>% 
  mutate(tag1=sapply(tag,function(x){
    p0 = str_split(x,"-",simplify = TRUE)[1,][1]
    p = str_split(p0,",",simplify = TRUE)[1,]
    start = str_split(x,"-",simplify = TRUE)[1,][2]
    stop = str_split(x,"-",simplify = TRUE)[1,][3]
    for(i in p){if(i>=start & i<=stop) return(1)}
    return(0)
  },simplify = TRUE)) %>% filter(tag1==1) %>% separate(col=tag,into=c('pos','start','stop'),sep = '-') %>% 
  unite("sequence_name",'chr','pos','MNV_id','ref','alt','MNV_rsid',sep = '_')
}


loss=left_join(ref,alt,by=c("sequence_name","motif_id",'motif_alt_id')) %>% 
  rename("ref_score"="score.x",
         "ref_p"="p-value.x",
         "alt_score"="score.y",
         "alt_p"="p-value.y",
         "ref_sequence"="matched_sequence.x",
         "alt_sequence"="matched_sequence.y",
  ) %>% 
  mutate(loss=sapply(alt_p,function(x){
    if(is.na(x))
      return("Loss")
    else
      return("None")
  } ,simplify = TRUE))

gain=right_join(ref,alt,by=c("sequence_name","motif_id",'motif_alt_id')) %>% 
  rename("ref_score"="score.x",
         "ref_p"="p-value.x",
         "alt_score"="score.y",
         "alt_p"="p-value.y",
         "ref_sequence"="matched_sequence.x",
         "alt_sequence"="matched_sequence.y",
  ) %>% 
  mutate(loss=sapply(ref_p,function(x){
    if(is.na(x))
      return("Gain")
    else
      return("None")
  } ,simplify = TRUE))

jmlist=read_table("/home/wangdy/MNVTFdb/Server/Services/Predict/resources/jaspar_motif_id",col_names = FALSE)
jmlist1=jmlist %>% mutate(X3 = toupper(X3))

if(nrow(loss)==0 && nrow(gain)>0){
  merge=gain %>% distinct() %>% inner_join(jmlist1,by=c("motif_id"="X2")) %>% 
    select(motif_id,motif_alt_id,X3,sequence_name,strand.x,ref_score,ref_p,start.x,stop.x,ref_sequence,strand.y,alt_score,alt_p,start.y,stop.y,alt_sequence,loss) 
} else if(nrow(loss)>0 && nrow(gain)==0){
  merge=loss %>% distinct() %>% inner_join(jmlist1,by=c("motif_id"="X2")) %>% 
    select(motif_id,motif_alt_id,X3,sequence_name,strand.x,ref_score,ref_p,start.x,stop.x,ref_sequence,strand.y,alt_score,alt_p,start.y,stop.y,alt_sequence,loss) 
}else {
  merge=gain %>% bind_rows(loss) %>% distinct() %>% inner_join(jmlist1,by=c("motif_id"="X2")) %>% 
    select(motif_id,motif_alt_id,X3,sequence_name,strand.x,ref_score,ref_p,start.x,stop.x,ref_sequence,strand.y,alt_score,alt_p,start.y,stop.y,alt_sequence,loss) 
  
}
  
write.table(merge,file = argv[8],quote = F,col.names = F,row.names = F,sep = "\t")
